FINITE DIFFERENCE WEIGHTS, SPECTRAL DIFFERENTIATION, 

AND SUPERCONVERGENCE 



Abstract. Let z\, Z2, ■ ■ . , Zn be a sequence of distinct grid points. A finite difference 
formula approximates the m-th derivative /( m )(0) as J2wkf (zk), with Wk being the 
weights. We derive an algorithm for finding the weights Wk which is an improvement of 
an algorithm of Fornberg (Mathematics of Computation, vol. 51 (1988), p. 699-706). 
This algorithm uses fewer arithmetic operations than that of Fornberg by a factor of 
4/(5m + 5) while being equally accurate. The algorithm that we derive computes finite 
difference weights accurately even when m, the order of the derivative, is as high as 
16. In addition, the algorithm generalizes easily to the efficient computation of spectral 
differentiation matrices. 

The order of accuracy of the finite difference formula for f( m \0) with grid points hzk, 
1 < k < N, is typically O (h N ~ m ) . However, the most commonly used finite difference 
formulas have an order of accuracy that is higher than the typical. For instance, the 
centered difference approximation (f(h) — 2/(0) + f(—h)) /h 2 to /"(0) has an order of 
accuracy equal to 2 not 1 . Even unsymmetric finite difference formulas can exhibit 
such superconvergence or boosted order of accuracy, as shown by the explicit algebraic 
condition that we derive. If the grid points are real, we prove a basic result stating 
that the order of accuracy can never be boosted by more than 1 . 



Since the beginning of the subject, finite difference methods have been widely used 
for the numerical solution of partial differential equations. Finite difference methods 
are easier to implement than finite element or spectral methods. For handling irregular 
domain geometry, finite difference methods are better than spectral methods but not as 
flexible as finite element discretizations. 

The basic problem in designing finite difference discretizations is to approximate 
/( m )(0), the m-th derivative of the function f(z) at z = 0, using function values at 
the grid points hzi, hz2, ■ ■ ■ , hzw- The grid points can be taken as zi, . . . , zn by setting 
the mesh parameter h = 1. We make the mesh parameter h explicit where necessary 
but suppress it otherwise. The finite difference formula can be given as either 
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1. Introduction 
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(1.2) 



f(m) (Q) 



W!,mf (hZi) + 



+ W N , m f (hz N ) 



h m 



NSF grant DMS-0715510 and SCREMS-1026317. 

1 



FINITE DIFFERENCE WEIGHTS 



2 



If we require (1.2) to have an error that is O (h N m ) for smooth /, the choice of the 
weights Wk, m , 1 < k < N, is unique (see Section 8). The grid points are always assumed 
to be distinct. 

Some finite difference formulas such as the centered difference approximations to 
f(0) and f»(Q)-(f(h)-f(-h))/2h and (f(h) - 2/(0) + f(-h)) /h\ respectively- 
occur very commonly and are part of the bread and butter of scientific computation. 
The most common finite difference formulas presuppose an evenly spaced grid. How- 
ever, evenly spaced grids are often inadequate. In applications, it is frequently necessary 
to make the grid finer within boundary layers or internal layers where the underlying 
phenomenon is characterized by rapid changes. In addition, evenly spaced grids do not 
lend themselves to adaptive mesh refinement. For these reasons, it is often necessary to 
use grids that are not evenly spaced. 

Fornberg [SJ QUI E] devised an algorithm for determining the weights Wi ym given the 
grid points z^. The results of this paper include two algorithms (see Sections 3 and 4) 
that improve Fornberg's. The numerical stability of algorithms to find finite difference 
weights can be subtle. Therefore we will begin this introduction by considering numerical 
stability. 

Numerical stability. All algorithms to compute finite difference weights come down to 
multiplying binomials of the form (z — Zk) and extracting coefficients from the product. 
The numerical stability of multiplying binomials, or equivalently of going from roots of a 
polynomial to its coefficients, has aspects that are not obvious at first sight. A dramatic 
example is the product (z — ou°) (z — u 1 ) . . . (z — u N ~ l ^j where u = exp(27ti/N). Math- 
ematically the answer is z N — 1. Numerically the error is as high as 10 15 for N = 128 in 
double precision arithmetic [4j. For numerical stability, the binomials must be ordered 
using the bit reversed ordering or the Leja ordering or some other scheme as shown by 
Calvetti and Reichel 0j. The roots must be ordered in such a way that the coefficients 
of intermediate products are not too large. 

Another point related to numerical stability comes up frequently. Suppose we want 
to multiply (z — a) into the polynomial (a® + ■ ■ • + auz M + ••■). The obvious way to 
form the coefficients of the product, which is denoted b + • • • + bj^z M + • • • , is to use 

bo = —aao 

(1.3) bj = -aaj + a,_i for 1 < j < M. 
A related problem is to assume a + a x z H = [z — a)^ 1 (b + b\z + • • • ) and find the 



a, in terms of the bj. The equations that comprise (1.3) can be easily inverted assuming 
a ^ 0: 

ao = — bo /a 

(1.4) a,j = {a,j-x — bj) / a. for j = 1,2,. 



Our point is that ( 1.3 ) appears to be safer numerically than ( 1.4 ). This point is discussed 



further in Section 6, but we note here that the system of equations (1.3) does not involve 



back substitution while the system of equations (1.4) involves back substitution. We 
codify this observation rule of thumb. 
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Rule of thumb: Algorithms that use triangular systems of recurrences with back substi- 



tution, for example systems such as (1.4), tend to be numerically unsafe. 

This rule of thumb will guide our derivation of a numerically stable algorithm for 
finding finite difference weights in Section 4. However, there are exceptions to it as we 
will see. 

Computation of finite difference weights. There exists a unique polynomial tt(z) of 
degree N — 1 which satisfies the interpolation conditions n(z k ) = fk for k = 1 . . . iV [BJ. 
The Lagrange form of this interpolating polynomial is given by 

TV 

(1.5) 7r(z) = ^w k n k (z)f k where n k (z) = J[(z - Zj) andw fe = l/n k (z k ). 

k=l j^k 

The finite difference weight w k<m is equal to the coefficient of z m in w k n k (z) times ml 
(see Section 8). The computation of the Lagrange weights w k takes 2N 2 arithmetic 
operations roughly half of which are multiplications and half are additions or subtractions 
(all operation counts are given to leading order only). 

In effect, Fornberg's algorithm [H] is to multiply the binomials (z — zj) using recursion 



of the form (1.3) to determine the coefficient of z m in the Lagrange cardinal function 
iu k 7v k (z k ). The algorithm is not presented in this way in [8j. Instead it is organized to 
yield the finite difference weights for partial lists of grid points z±, . . . ,z k with k increasing 
from 1 to N. Fornberg's algorithm requires 5N 2 /2 + 5MN 2 /2 arithmetic operations to 
determine the weights w kjTn , 1 < k < N, if m = M. In the operation count, the 
coefficient of N 2 is proportional to M because each Lagrange cardinal function w k 7i k (z) 
is treated independently. Since the order of the derivative goes up to M only, coefficients 
beyond the z M term are not needed and are not computed. 

The algorithm for finding finite difference weights presented in Section 3 uses the 
modified Lagrange formula [2]: 

(1.6) tt(z) = {z- Zl )...(z- z n ) f-^-A + • • • + -^-fjr) • 

\Z — Zi z — Z N J 

Here the idea is to begin by determining the Lagrange weights w k and the coefficients (up 
to the z M term) of the polynomial n*{z) = Y\ k=1 (z — z k ) . The coefficients are determined 



using a recursion of the form (1.3) repeatedly. Since w k Tt k (z) = w k n*(z)/(z — z k ) , we 



may then use a recursion of the form (1.4) to determine the finite difference weights. 
This algorithm uses 2N 2 + 6MN operations. The 2N 2 term is the expense of computing 
the Lagrange weights w k . 



Since the method based on the modified Lagrange formula uses (1.4), it is numerically 
unsafe according to our rule of thumb. Indeed, for large M it is not numerically stable. 
The computations of Section 7 suggest that the method based on the modified Lagrange 
formula is a good choice for M < 4 but not for larger M. 

In Section 4, we derive another algorithm, one based on the following partial products: 



k N 

h(z) = Y[(z - zj) and r k (z) = J[(z - Zj). 
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By convention, l = r^ + i = 1. For k = 1, . . . , N, coefficients of the se partial products 
are computed up to the z M term using recursions of the form (1.3) repeatedly. Since 
7Tfc = Ik-iTk+i, the finite difference weights Wk, m for m = 0, . . . , M, are obtained by 
convolving the coefficients of 4_i and r& + i followed by a multiplication by ml and the 
Lagrange weight w^- This algorithm uses 2N 2 + 6NM + NM 2 arithmetic operation. 

Even though the method based on partial products has an additional expense of NM 2 
operations, we recommend it for all uses. Because it completely avoids back substitution, 
it has good numerical stability. If FFTs are used for convolution, the NM 2 term can 
be replaced by O (N M log M) . Instances where M is so large that the use of FFTs is 
advantageous are unlikely to occur in practice. 

Spectral differentiation. Beginning with ( 1.1 ), our discussion of finite difference weights 
has assumed z = to be the point of differentiation. Given grid points Zi, . . . , Zjy, the 
spectral differentiation matrix of order M is an N x N matrix. If the (i, j)-th entry is 
denoted Uij, then f( M \zi) « Y^j^ijfizj). Thus coij is the finite difference weight at Zj 
if the point of differentiation is z = z^. The point of differentiation can be shifted to 
by replacing the grid points zi, . . . , z^ by z\ — z i: . . . , z^ — z± . 

Applying Fornberg's method row by row would cost O (N 3 ) arithmetic operations. 
Welfert [18J modified Fornberg's recurrences and obtained a method that computes the 
spectral differentiation matrix of order M using only O (N 2 M) operations. 

The way to modify the first of our two methods (Section 3) so as to compute spectral 
differentiation matrices in O (N 2 M) arithmetic operations is almost obvious. We simply 
have to note that the Lagrange weights Wk defined by ( 1.5[ ) do not change at all when 
the entire grid is shifted. Therefore Lagrange weights are the same for every row of 
the spectral differentiation matrix and need to be computed just once using O (N 2 ) 
operations. The rest of the computation is repeated for every row, with an appropriately 
shifted grid, costing O (N 2 M) operations to determine the entire spectral differentiation 
matrix. 

If the method based on partial products (Section 4) is used to determine the coeffi- 
cients of the Lagrange cardinal functions Wk^k( z )i the cost of computing the spectral 
differentiation matrix is O (N 2 M 2 ). Although Welfert's method and the method of Sec- 
tion 3 compute spectral differentiation matrices with a lower asymptotic cost, they are 
less accurate than the method based on partial products. They should not be used for 
M > 4. Welfert [TS] stated that the problem of round-off errors becomes "very im- 
portant" for M > 6 and that his tables did not use a large enough M to expose the 
problem. In contrast, the method based on partial products computes every entry of 
the 512 x 512 Chebyshev differentiation matrix for the 16-th derivative with 9 or more 
digits of accuracy. 

Superconvergence or boosted order of accuracy. If the number of grid points is iV and 
the order of the derivative is m, the finite difference weights are unique if the difference 
formula is required to be O (h N ~ m ^j (see Section 8). For certain grids, these unique 

weights imply an error of O (ji N ~ m+1 ^J, which is of higher order than what is typical for 
iV grid points and the m-th derivative. We term this as superconvergence or boosted 
order of accuracy. In section 8, we gi ve e xplicit conditions for boosted order of accuracy. 
The finite difference approximation (1.2) to p m ^(0) has an order of accuracy boosted by 
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1 if and only if 

where Sk is the elementary symmetric function 



E 



z i\... z i k ■ 



l<h<-<i k <N 

If the grid points are real, we prove that the order of accuracy cannot be boosted by 
more than 1. 

For the special case m — 2, the finite difference approximation to the second derivative 
at z = using three grid points has an order of accuracy equal to 2, which is a boost of 1, 
if and only if the grid points satisfy z\ + z 2 + z 3 = 0. Evidently, this condition is satisfied 
by the grid points —1,0,1 used by the centered difference formula. An unsymmetric 
choice of grid points such as —3, 1, 2 also boosts the order of accuracy by 1. However, 
no choice of Zi, Z2, and z^ on the real line can boost the order of accuracy by more than 
1. 

With four grid points and m = 2, the condition for a boost in the order of accuracy is 

Z 1 Z 2 + Z\Zz + Z\ Z A + Z 2 Z3 + Z 2 Z 4 + Z Z Z± = 0. 

No choice of z±, z 2 , z 3 , and z 3 on the real line can boost the order of accuracy of the 
finite difference approximation to /"(0) by more than 1. The maximum possible order 
of accuracy is 3. In this symmetric choice of grid points such as —2,-1,1,2 

does not boost the order of accuracy. Unsymmetric grid points that boost the order of 
accuracy to 3 can be found easily. For example, the order of accuracy is 3 for the grid 
points -2/3,0,1,2. 

These results about super convergence or boosted order of accuracy of finite difference 
formulas are quite basic. It is natural to suspect that they may have been discovered a 
long time ago. However, the results are neither stated nor proved in any source that we 
know of. 

If the grid points Z{ are allowed to be complex, the order of accuracy can be boosted 
further but not by more than m. The order of accuracy is boosted by k with 1 < k < m 
if and only if 

Sn-tu = dN-m+l = • • • = SjV-m+fc-1 = 0. 

An algorithm to detect the order of accuracy and compute the error constant of the 



finite difference formula (1.2) is given in section 8. 



2. From roots to coefficients 

Given a±, . . . , a^, the problem is to determine the coefficients of a polynomial of degree 
iV whose roots are cnx, • • • , (%n- The polynomial is evidently given by Y\^ =1 (z — a^). If 
the product II/Li^ — a k) is given by c + C\Z + ■ ■ ■ + z n , then the coefficients c' , c[, . . . 
of the product Ilfcii(<2 — a k) are formed using c' = — c a„ + i and c' m = —c m a n+ i + c m _i 



for m = 1, 2, ... as in (1.3). 



For accurate evaluation of the coefficients, Calvetti and Reichel [1] have demonstrated 
the need to order the roots carefully. The roots must be ordered in such a way 
that the coefficients of the partial products I\k=i( z ~ a k), n — 1, . . . , N — 1, that occur in 
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intermediate stages are not much bigger than the coefficients of the complete product. If 
aj- = where oj = exp(27ri/iV), the natural ordering oj°, oj 1 , . . . , oj n_1 is numerically 

unsound. If iV is large the first roots in this sequence are close to 1 leading to partial 
products which resemble (z — l) n and have coefficients that are of the order of the 
binomial coefficients. In contrast, the complete product is simply z — 1. 

This matter of ordering the roots carefully is equivalent to choosing a good order of 
grid points when determining finite difference weights. Good ordering of grid points may 
improve accuracy but is not as important as it is in the general problem of determining 
coefficients from roots. When determining finite difference weights for a derivative of 
order M, we need coefficients of terms 1, z, . . . , z M but no higher. The most dramatic 
numerical stabilities in determining coefficients occur near the middle of the polynomial, 
but M, which is the order of differentiation, will not be large in the determination of 
finite difference weights. 

Suppose that the «j are all nonzero and that T\k=i( z ~ a k) = (% + ••• + cmz m + 
O (z M+1 ^j . Another algorithm to compute Cq,...,Cm is obtained as follows. Let 

N 

v r = E«r r 

k=l 

£ r = Y, (a h . . . a^y 1 ■ 

l<ii<—<i r <N 

By the Newton identities 

£\ = Vi 
2£ 2 = £ 1 V 1 -P 2 
3£ 3 = £ 2 V 1 -£ 1 V 2 + V 3 



and so on. By convention, £ = 1. The algorithm begins by computing the power 
sums Pi,..., Vm directly and uses the Newton identities to compute the elementary 
symmetric functions £ r , < r < M. The coefficients are obtained using 



N 

Cr = (-l) N ^-%l[a k . 

k=l 



This algorithm does not really presuppose an ordering of the ak and the computation 
of the power sums V r is backward stable, and especially so if compensated summation 
is used [12] . If this method is used to compute the product n^Li ( z ~ ujk1 ^ where oj is 
as before, it finds the coefficients of the product with excellent accuracy. But in general 
this method is inferior to the repeated use of (1.3) after choosing a good ordering of the 
roots By way of a partial explanation, we note that the Newton identities have a 
triangular structure with back substitution, which is deemed to be possibly unsound by 
the rule of thumb stated in the introduction. 
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3. Finite difference weights using the modified lagrange formula 

, zn with fx, . . . , /at being the function values at the grid 



Let the grid points be Z\ 
points. Define 

(3.1) 



n k {z) 



N 

n< 



Z — Za 



Then the Lagrange interpolant shown in (1.5) is tt(z) = J2k=i w k 7T k( z )fk- The weights 
w k equal l/u k {z k ). Our objective is to derive formulas for d m 7r(z)/dz m at z = for 
m = 1,...,M. The m = case is regular Lagrange interpolation. The weights w k 
will be assumed to be known. The formulas for d Tn 7c(z)/dz m at z = will be linear 
combinations of f k with weights. We assume 1 < M < N — 1. 
If the coefficient of z m in iTk(z) is denoted by Ck, m , we have 



d rn 7i(z) 
dz m 



2 = 



N 

m\ Ck,mW k fk- 
k=l 



The finite difference weights are then given by 
(3.2) w ktTn = m\w k c k>m . 



Once the Ck, m are known, the weights w k ,m, are computed using (3.2) for k 
and m = 1, ... , M. 

Let 7T*(z) denote the polynomial IIfe=i( 2; — z k 

). Let 



1,...,N 



71 



N 

[z) = E C k z k . 

k =0 



Notice that tt*{z) occurs as a factor in front of the modified Lagrange formula (1.6). Our 
method for calculating w k ^ m begins by calculating C , Ci, . . . , Cm+i- The c k . m are deter- 
mined using C , . . . , Cm+i and then the w ktTn are determined using (3.2). The Lagrange 



weights figure in this last step. It is in this sense that the method for determining finite 
difference weights described in this section uses the modified Lagrange formula. 



We start by setting Co = 1 and C\ 
following update is performed: 



a 



= Cm+i = 0. For each j = 1, 2, . . . , iV , the 

-ZjCq 

-ZjCi + Cq 



M 



Ci 



M+l 



-ZjCM + Cm-1 



-ZiC 



■jls M+l 



+ c 



M 



followed by C = C , ... , C M +i = C' M+1 . 
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To obtain the Ck >m , use (z — Zk) 7Tk(z) = 7r*(z) to get 

—ZkCkfl — Co 



If Zk 7^ 0, we have 



-%k c k,M + Cfc,A/-l — Cm 
-Zk c k,M+l + c k,M — Cm+1- 

c k,0 — —Cs/Zk 

Ck,i = {ck,o-Cx)/zk 



(3.3) Ck.M — {c-k.M-i — Cm) I ' Zk- 

If Zk = 0, we have Ck )m = C m+ i for m — 0, 1, . . . , M. We can now use (3.2) to find the 
finite difference weights. The complete algorithm is exhibited as Algorithm [1} 

The operation count for Algorithm [T] is as follows. The operation counts are given to 
leading order only. 

• The function LagrangeWeights() invoked on line 34 computes the Wk using 2N 2 
operations. More precisely, the operations are N(N — 1) subtractions, N(N — 2) 
multiplications, and iV divisions. The number of subtractions can be halved 
using additional storage. 

• The function FindCQ invoked on line 36 computes Co, . . . , Cm+x- The num- 
ber of operations used is 2MN. More precisely, the operations are N(M + 2) 
multiplications and N(M + 1) additions. 

• The function findckmQ invoked N times on line 39 computes Cfc m . The num- 
ber of operations used is 2MN. More precisely, the operations are N(M + 1) 
multiplications, NM additions, and one division. 

• The function Find Weights () invoked iV times on line 40 computes the finite 
difference weights Wk, m - The number of operations used is 2MN. This function 



implements Wk, m = f^-WkCk,m-, which is (3.2), using two multiplications for each 
Wk, m - Even if the factorials ml are precomputed and stored, we need the same 
number of multiplications for each Wk, m 
The total number of floating point operations is 2N 2 + 6MN. 

The operation count of Fornberg's method is 5N 2 /2 + 5MN 2 /2 — 5M 3 /6 assuming 
M -C N. Our algorithm differs from that of Fornberg [H] in two major respects. Firstly, 
Fornberg does not compute the Lagrange weights Wk explicitly as we do. Secondly, 
we form the coefficients of tt*{z) and use that to recover the coefficients of n k (z) for 
k = 1, . . . , N. Fornberg's method is laid out quite differently from ours, but in effect it 
treats each n k (z) separately. 

Because Fornberg's method builds up the finite difference weights for the grid zi, . . . ,z N 
using the finite difference weights of the partial grids Z\ , . . . , Z\. , with k increasing from 
1 to N, it is forced to use O (N 2 ) divisions. Algorithm [l] uses only 2N divisions. These 
occur in the computation of the Lagrange weights (line 4) and in determining Ck im (line 
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Algorithm 1 Finite difference weights using the modified Lagrange formula 



1: function LagrangeWeights(zi, . . . , Z N ,Wi, . . . , w N ) 
2: for % = 1,2,..., TV do 

Wi = l\j( z i — z j) over J 1 — 1) • • • , N but j ^ i. 



3 
4: 
5 
G 
7: 
8 
9 
10 
11 
12 
13 
14 
15 
1G 
17 
18 
19 
20 
21 
22 
23 
24 
25 
2G 
27 
28 
29 
30 
31 
32 
33 
34 
35 
3G 
37 
38 
39 
40 
41 
42 



Wi = 1/Wi 

end for 
end function 

function FindC(zi, . . . , zn, C , . . . , Cm+i) 
Temporaries: to, • • • , £m+i 
C = 1 and d = for 1 < % < M + 1 
for j — 1, . . . , N do 

to — —ZjCo 

ti = Ci_i - zjCi for i = 1, 2, . . . , M + 1 
Ci = ti for i = 0,l,...,M + l 
end for 
end function 

function FlNDCKM(z fc , C , Ci, . . . , C M +i, c fc>0 , • • • , c fcjM ) 
if £ fc == then 

c fc ,m = C m+ i for m = 0, . . . , M 
else 

C = l/Zk 

c k ,m = C(cfc,m-i - C fc ) for m = 1, . . . , M 
end if 
end function 

function FlNDWElGHTs(c fci0 , • • • , Cfc,M, w k , w kfi , w ktM ) 
f = w k 

for m = 0, 1, .. . ,M do 

^fc,m fc k ,m 

f = (m+l)f 
end for 
end function 

function FindAllWeights(zi, ...,z n , w k , m for k = I, . . . , N and m = 0, . . . , M) 

Temporaries: Wi, . . . , Wjy 

Lagrange Weights(^i, . . . , z N ,w!, w N ) 

Temporaries: Co, . . . , Cm+i 

FindC(2;i, . . . , z N , C , . . . , Cm+i) 

Temporaries: c k)m for k = 1, . . . , iV and m — 0, . . . , M 

for k = 1, . . . , N do 

findckm(^, C , Ci, . . . , C M +i, c kfi , . . . , c kjM ) 
FindWeights(c m , • • • , c kM , w k , w kfi , w KM ) 

end for 
end function 
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20). Similarly, Algorithm [2j which is derived in the next section, uses only N divisions. 
On current processors, division is more expensive than multiplications or additions. For 
example, in the Intel Nehalem microarchitecture, the latency of division is three to six 
times that of multiplication. While multiplication instructions can be dispatched to 
ports in successive clock cycles, the dispatch of division instructions must be separated 
by five or so clock cycles. 

Many, if not most, of the numerical analysis textbooks recommend the Newton form 
for polynomial interpolation. The weights of the Newton form can be computed with 
O {N 2 ) arithmetic operations using the divided differences table and updated using 
O(N) operations if a new grid point z^ + i is added. In addition, the Newton form can 
be evaluated at a point z using O(N) operations. It has been well known that the 
weights of the Lagrange form can be computed with O (N 2 ) operations, but there was 
much less clarity about the evaluation of the Lagrange form and the cost of updating 
the weights when a new grid point z^ + i is added. In an engaging paper, Berrut and 



Trefethen j2] pointed out that an examination of (1.6) and the formula (1.5) for the 
Lagrange weights Wk clarifies the updating and evaluation of the Lagrange form to be 
as efficient as in the Newton case. For the use of the Lagrange form for finding roots of 
functions, see Corless and Watt [5]. 

4. Finite difference weights using partial products 

As we will see in Section 7, Algorithm [T] is accurate enough if M < 4, but it should 
not be used if the order of the derivative is higher than 4. The problem is the use of 



(3.3 ) by the function findckm ( ) to determine c k>m . This step involves back substitution. 
We will now derive an algorithm that completely avoids back substitution. 

Let lk(z) = Ylj=i{ z — z j) and r k (z) = Y\f =k (z — Zj). Denote the coefficients of 
1, z, . . . , z M in l k (z) and r k (z) by L kj0 , ■ ■ ■ , L k]M and R kt0 , . . . , R k ,M, respectively. The 
coefficients Lk, m are computed in the order k = 1,2, . . . , N. The coefficients Rk, m are 
computed in the reverse order, which is k = N, N — 1, . . . , 1. It is evident that ir k (z), 
which is defined by ( |1.5 ) or (3.1), is equal to l k -i(z)r k+ i(z). Therefore the coefficient 



Ck : m of z m in ir k (z) can be obtained using 

m 

Cfc,m ^ ; Lk—l t m—sRk+l,S' 

The finite difference weight w k ,m is obtained as m\wkCk, m , where w k is the Lagrange 
weight at Zk- 

Because this method stores L k ^ m , Rk t m, and other intermediate quantities, it is more 
convenient to implement it as a class than in purely functional form. The public members 
of the C++ class are shown below: 

class FDWeights{ 

FDWeights( const double *zz , int NN, int MM) ; 

-FDWeights ( ) ; 

void setzO (double zO ) ; //derivative at zO 

void setzk(int kO ) ; //derivative at kO—th grid point 

double operator ()( int m, int k){ //weight for mth derv at kth grid point 
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assert ((0<=m)&&(nK^&&(0<=k)&&(k<N) ) ; 
return fdw [k*(Mfl)+m] ; 

} 

double operator ()( int k){//weight for **Mth** derv at kth grid point 
assert ((0<=k)&&(k<N) ) ; 
return fdw [k*(Mfl)+M] ; 

} 

}; 

The fdw array stores the finite difference weights so that fdw [k* (M+l) +m] is equal 
to Wk,m- The implementation of the member functions is displayed as Algorithm [2] . 

Algorithm [2] invokes functions defined as a part of Algorithm [I] on lines 8 and 35. The 
total expense is 2N 2 + 6NM + NM 2 arithmetic operations. 



5. Spectral differentiation matrices 

In Algorithm [2j the Lagrange weights Wk are computed in the class constructor 
FDWeights: :FDWeights () . If we want the finite difference weights for derivatives eval- 
uated at z = £, we need to invoke the member function FDWeights: :setzO () with the 
argument (. If we want the finite difference weights for the derivatives evaluated at 
the fc-th grid point Zk, we need to invoke the member function FDWeights: :setzk(). 
The Lagrange weights are not re-computed when either of these member functions is 
used to set the point at which derivatives are taken. This has implications for spectral 
differentiation. 

Suppose the grid to be zi,...,z N as usual. In the N x N spectral differentiation 
matrix of order M, the (i,j)-ih entry Uij is equal to the finite difference weight at Zj 
when the derivative of order M is taken at z — Zi. The spectral differentiation matrix 
is computed in the following steps. 

l: FDWeights fd(si, . . . , z N ,M) 
2: for i = 1, . . . , N do 
3: fd.setzk(z) 

4: Wij=fdtf) for j = 1,...,N 
5: end for 

In the class definition of the previous section, the function call operator has been over- 
loaded so that f d ( j ) returns the finite difference weight at Zj for the M-th derivative. 

Since the Lagrange weights are computed just once, the cost of computing the spectral 
differentiation matrix of order M is 2N 2 + QN 2 M + N 2 M 2 arithmetic operations. 



6. Discussion of numerical stability 
Suppose that the sequence a , a±, . . . and the sequence b , bi, . . . are related by 

(a + a x z + a 2 z 2 H ) = (z - a) -1 (b + b x z + b 2 z 2 H j . 
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Algorithm 2 Finite difference weights using partial products. 

1: function FDWEIGHTS::FDWEIGHTS(zi, . . . , z N ,M) 
2: Variables internal to class 

3: (1): Order of derivative M (initialized from argument list) 
4: (2): Grid points zi, . . . ,Zn (initialized from argument list) 
5: (3): Lagrange weights wi, . . . , wn 

6: (4): Partial product coefficients: L km and Rk, m for 0<k<N + l,0<m<M 
7: (5): Finite difference weights: w k ^ m for 1 < k < N and < m < M 
8: Lagrange Weights^, . . . , z N ,wi, . . . , w N ) 
9: setzO(O) (derivatives at z = by default) 

10: end function 

11: function FDWEIGHTS::-FDWEIGHTS 
12: Deallocate all internal variables 

13: end function 

14: function MULTBlNOM(a , . . . , a,M,bo, • • • , &m,C) 
15: b = -( a 

16: b rn = -(a m + a m _i for k = 1, . . . , M 

17: end function 

18: function cONVOLVE(a , . . . , a M ,b , . . . , b M ,c , . . . , c M ) 

19: c m = a m b + a m _i6i H h a b m for m = 0, . . . , M 

20: end function 

21: function FDWEIGHTS::SETZO(C) 

22: Temporaries: (i, . . . ,(n 

23: ( k = z k -( for k = l,...,N 

24: L m = 1 for m = and L fe m = for m — 1, . . . , M 
25: for k = 1,...,N do 

26: MULTBINOM(L fc _ li0 , . . . , L k _ 1>Ml L kfi) L k;M , ( k ) 

27: end for 

28: RN+i,m = 1 for m = and Rn+i,™, = for m = 1, . . . , M. 
29: for k = N, N - 1, . . . , 1 do 

30: MULTBINOM(i? fe+ i )0 , . . . , R k+1>M ,Rkfi, • • • , Rk,M,Ck) 

31: end for 

32: for fc = 1, . . . , N do 

33: Temporaries: c k>m 

34: CONVOLVE(L fc _ li0 , . . . , L k _ 1;M ,R k+lj0 , . . . , R k +i,M,C k fi, • • • , c A:,m) 

35: FlNDWEIGHTS(c M , . . . , C fcjM ,ti; fc ,W M , • • • , Wfc,m) 

36: end for 

37: end function 

38: function SETK(k) 
39: SETZO(Zfc) 

40: end function 
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Then the relationship between the sequences can be expressed in matrix notation as 
(6.1) 

/ -a \ / a \ ( b \ ( ~\ i \ f b ° \ f a ° \ 

ai 



1 -a J 



\ b M J 



a 
a i 



IS 



\ a M + l ■ a 2 a / \ / 

The matrices that occur here will be denoted by T and T -1 , respectively. To calculate 
the a,j given the bj, the method used ordinarily is back substitution — a = —b /a and 
cij = (cij-i — bj)/a for j > 1 — as in (1.4). However, in Section 1, we stated as a rule of 
thumb that an intermediate step which uses back substitution of that type is likely to be 
numerically unsafe. In particular, that rule of thumb leads us to expect that Algorithm 
[TJ which uses back substitution (to compute Ck, m ), is possibly inferior to Algorithm [2j 
which does not use back substitution. In the next section, we will show that that surmise 
is indeed true. Here we will discuss the rule of thumb. 

For convenience, we denote the two vectors that occur in (6.1) as a and b. If a 
calculated using back substitution, the following norm-wise bound on errors applies 
Chapter 7]: 

||a-a|| 2e«(T) 
a - 1 -€k(T)' 

Here a is the computed vector, e is a small multiple of the machine epsilon, and k(T) = 
||T|| | (T" 1 1 1 is the condition number. The norm can be any matrix norm such as the 
2-norm or the oo-norm. 

By inspecting T and T~ l displayed in (6.1), it is evident that k(T) increases expo- 
nentially with M if | a | < 1. The norm-wise bound suggests that the norm- wise error in 
a increases exponentially with M. 

If | a | > 1, the matrix T is evidently well-conditioned. So it may appear as if the 
problem can be cured by recasting it to make \a\ > 1. Such a recasting is easy to 
accomplish. If we write (z — a) = d(z/d — a/d) for some d < \a\ and expand the series 
in powers of {z/d), then a/d will replace a in the triangular systems and T will therefore 
be well-conditioned. Such scaling does not improve accuracy, however. It is true that the 
norm-wise errors in the computed vector will be small, but the entries of the computed 
vector will be poorly scaled. When the computed vector is multiplied by inverse powers 
of d to recover entries of a, the relative errors in entries such as au can be quite large. 

Since T is triangular, a component-wise bound of the following type applies [T21 Chap- 
ter 81: 



\ a M J 



a 




a 


|oo 




a 


| oo 



< 



cond(T, a)7„ 



-, cond(T, x) 





T -l 


1 \ T \ 


X 








x\ 


\ oo 



cond(T, x) 



\T\ IT" 



1 — cond(T)7„ 

Here 7„ is approximately 2n times the machine epsilon and |T| is T with its entries 
replaced by their absolute values. This component-wise bound does not help much. If 
| a | < 1, the condition numbers that occur in this bound again increase exponentially 
with M. On the other hand, if the power series are expanded in powers of z/d for some 
d < | a |, the condition numbers become mild but the bound on the oo-norm relative 
error becomes useless when the computed vector is scaled by inverse powers of d. 
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In further support of the rule of thumb, we mention that triangular matrices are 
typically ill-conditioned and their condition number increases exponentially with the 
dimension of the matrix [TB]. For example, if all entries of a triangular matrix are 
independent normal variables with mean 0, the condition number increases at the rate 
2 M , M being the dimension of the matrix. 

An example where the accuracy does not deteriorate rapidly with M in spite of the use 
of back substitution in an intermediate step occurs in barycentric Hermite interpolation 

Triangular systems with back substitution come up in a natural way if we want to 
determine the coefficients bj such that 

b + b x z + b 2 z 2 + ■■ 1 



a + a\z + a 2 z 2 H ' 

However, the discussion here suggests that such a method will be inaccurate. The 
fast Fourier transform (FFT) and contour integrals may have a role in the numerically 
accurate invertion of series — see [3] for related ideas. 

7. Numerical Examples 

For simple choices of grid points, such as Zk = 0, ±1, ±2, ±3, ±4, Algorithms [T] and [2j 
which are respectively based on the modified Lagrange formula (Section 3) and partial 
products (Section 4), as well as Fornberg's method find the finite difference weights with 
errors that are very close to machine precision. To compare the different methods, we 
must turn to more complicated examples. 

The Chebyshev points are defined by 

z k = cos ((jfe - 1)tt/(JV - 1)) = sin (vr(iV - 2k + 1)/(JV - 1)) 

for k = 1, . . . , N. We will look at the relative errors in the spectral differentiation 
matrix of order M for M = 2, 4, 8, 16 and N = 32, 64, 128, 256, 512. When Algorithm [2] 
is employed, the spectral differentiation matrix is computed as described in Section 5. 

The Chebyshev points are distributed over the interval [—1, 1]. The logarithmic ca- 
pacity of an interval is one quarter its length, which in this case is 1/2. Therefore the 
Lagrange weights Wk will be approximately of the order 1/2 N . To prevent the possibility 
of underflow for large N, the Chebyshev points are scaled to 2zk and the resulting finite 
difference weights for the M-th derivative are multiplied by 2~ M . 

For reasons described in Section 2, the Chebyshev points are reordered. The reorder- 
ing we use is bit reversal. With N being a power of 2 in our examples, the binary 
representation of k (here k is assumed to run from 1 to N — 1) can be reversed to map 
it to a new position. The permutation induced by bit reversal is its own inverse, which 
simplifies implementation. The reordering of the grid points has the additional effect of 
making underflows less likely |14j . For a discussion of various orderings of Chebyshev 
points, see jl] . The figures and plots here are given with the usual ordering of Chebyshev 
points. 



Before turning to Figures 7.1 and 7.2, which compare the numerical errors in different 
methods, we make an important point. Even though the number of digits of precision 
lost in the 32 x 32 differentiation matrix of order M = 8 may be just 3, the errors in 
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(a) Algorithm [l] (b) Algorithm |2] (c) Fornberg's algor 

Figure 7. 1. Errors in the entries of the 32 x 32 Chebyshev differentiaton 
matrix of order M = 8. The vertical axis labeled d shows the number of 
digits of precision lost due to rounding errors. 



an 8-th derivative evaluated using that matrix will be much higher. Some entries of the 
N x N Chebyshev differentiation matrix of order M are O (^N 2M ^j. Very large entries 
occur in the differentiation matrix and in exact arithmetic an accurate derivative will 
be produced after delicate cancelations during matrix-vector multiplication. In finite 
precision arithmetic, the largeness of the entries implies that even rounding errors in the 
entries that are of the order of machine epsilon are sufficient to cause explosive errors in 
numerically computed derivatives. The 512 x 512 Chebyshev differentiation matrix of 
order 16 is useless even if every entry is computed with the maximum possible 16 digits 
in double precision arithmetic. 

There are tricks for improving the accuracy of computed derivatives or of entries 
of the Chebyshev differentiation matrix [T7j. The so-called negative sum trick can 
be interpreted as a barycentric formula for a derivative [U ITS] . This trick is limited 
to derivatives of the first order, and even when M — 1, it does not help when the 
differentiation matrix is inverted in some form. 

Our purpose here is to assess the accuracy with which the finite difference weights are 
computed and we will stick to that purpose. From Figure 7.1, we see that Algorithm [TJ 



which is based on the modified Lagrange formula, loses 7 digits for N = 32 and M 
while the other two methods lose only 3 digits. There is a kind of flip symmetry in the 
errors shown in each of the plots of that figure. 

Figure 1_J2_ gives a more extensive report of errors. All the errors were estimated using 
50 digit arithmetic in MAPLE The errors were validated using 60 digit arithmetic. From 
the figure, we see that the algorithm based on partial products (A2 in the legends of 
the figure) is as accurate as Fornberg's method in spite of using many fewer arithmetic 
operations. From the four plots of Figure 7.2, a surprise is that the errors are smaller 
for M = 16 than for M = 4 or M = 8. Why is that the case? We are not certain of the 
answer. 
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10* 10 2 

N N 

(c) M = 8 (d) M = 16 



Figure 7.2. Variation of the maximum relative error over the N 2 entries 
of the NxN Chebyshev differentiation matrix. The order of differentiation 
is M. In the legends, Al and A2 stand for Algorithms [T] and [2j respectively. 

The errors in the algorithm based on the modified Lagrange formula (Al in the leg- 
ends) is already noticeably larger for M — 8. For M = 16, they are quite bad and of 
the order of 10 5 (omitted from Figure |7.2p ). Such a rapid deterioration in error with 
increase in M appears to validate the exponential instability phenomenon hinted at in 
the previous section. 

The main finding of this section is that Algorithm[2j which is based on partial products, 
is as accurate as Fornberg's method even though it uses fewer arithmetic operations. 



8. SUPERCONVERGENCE OR BOOSTED ORDER OF ACCURACY 

Let zi, . . .Zn be distinct grid points. Let 



.r 



f(m) (Q) 



Wl,mf (hz{) j h W NiTn f (hz N ) 

h m 



be an approximation to the m-th derivative at 0. We begin by looking at the order of 



accuracy of this approximation. Here (1.2) is shown again as (8.1) for convenience. The 
order of the derivative m is assumed to satisfy m < N — 1. The case m = corresponds 
to interpolation. The allowed values of m are from the set {1, 2, . . . , N — 1}. 
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Lemma 1. The finite difference formula (8.1) has an error of O (h 

N 

and w k,mXk = 



N-m 



if and only if 



N 



k=l 



ml 



k=i 



for n G {0, 1, . 
differentiable. 



N — 1} — {m}. The function f is assumed to be N times continuously 



Proof. Assume that the weights Wk, m satisfy the conditions given in the lemma. The 

function f(z) can be expanded using Taylor series as f(0)+f'(0)z-\ \-f ( * N ~ 1 \0)z N ~ 1 / (N- 

1)! + z N g(z), where g{z)is a continuous function. In particular, g(z) is continuous at 
z = 0. If the Taylor expansion is substituted into the right hand side of (8.1) and the 



conditions satisfied by the weights are used, we get the following expression: 
/M( ) + h N ~ m (w 1>m z?g (h Zl ) + ■■■ + w N , m z%g (hz N )) . 

The coefficient of h N ~ m is bounded in the limit h — > 0, and therefore the error is 
0{h N ~ m ). 

The necessity of the conditions on the weights Wk, m is deduced by applying the finite 
difference formula (8.1) to / = 1,2, . 



JV-l 



□ 



The conditions on the weights in Lemma [3] correspond to the following matrix system. 



(8.2) 



1 



1 

z 2 



. y N-l 7 JV-1 



Z 



1 

Z N 

N-l 
N 



\ 



( Wi,m 



\ 



m\e r 



\ U)N,m J 



where e m is the unit vector with its m-th entry equal to 1. The matrix here is the 
transpose of the well-known Gram or Vandermonde matrix. 

Newton and Lagrange interpolation are techniques for solving Vandermonde systems. 
Newton interpolation is equivalent to an LU decomposition of the Gram or Vander- 
monde matrix [B]. Partly because the matrix in (8.2) is the transpose of the Gram or 
Vandermonde matrix, the interpolation techniques are not directly applicable. 

The Gram or Vandermonde determinant equals ni<i<j<jv ( z j ~ z i) an d is therefore 
nonsingular [UJ. Thus we have the following theorem. 

Theorem 2. There exists a unique choice of weights Wk. m , k = 1, . . . ,N, such that the 
finite difference formula (8.1) has error O (h N ~ m 



This theorem is trivial and generally known. However, its clear formulation is essential 
for developments that will follow. Our main interest is in boosted order of accuracy. 



Lemma 3. The finite difference formula (1.2) has boosted order of accuracy with an 
error of O (h N ~ m+b ) , where b is a positive integer, if and only if the weights Wk jm satisfy 



N-i+p 



i N-1+/3 n 

W\,m z l H h WN,m z N = 

for /3 — 1, . . . , b in addition to the conditions of Lemma [IJ 
Proof. Similar to the proof of Lemma [TJ 



□ 
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To derive conditions for boosted order of accuracy that do not involve the weights, 
we introduce the following notation. By 



(8.3) 



det (zi, z 2 ... z N ; m, n 2 , • • • , n N ) 



we denote the determinant of the NxN matrix whose (i,j)-th entry is z™\ The transpose 

of the grid points, which 

Zn; 0, ... N — 1) in this notation 
Theorem 4. Let Wk m , k 



of the Vandermonde or Gram determinant of the grid points, which occurs in (8.2), is 
det(zi, 



1, . . . , N, be the unique solution of (8.2) so that the finite 
difference formula (8.1) has an order of accuracy that is at least N — m. The order of 
accuracy is boosted by b, where b is a positive integer, if and only if 

det (zi, ...,z N ][Q,l,...,N-l,N-l + p]-m) = 

for /3 — 1, . . . , b. Here [0, 1, . . . , N — 1, N — 1 + 0\ — m denotes the sequence 0,1, ... N — 
1, N — 1 + /3 with m deleted. 

Proof. First, assume the weights Wk, m and the grid points Zf. to be real. The condition 
of Lemma |3| requires that the row vector W m = [wi m , . . . , WN,m\ be orthogonal to 



•4) 



z l i 



Z 



N 



By (8.2), W m is orthogonal to every row of the Gram matrix except the m-th row. Since 
the Gram matrix is non-singular, the rows of that matrix are a linearly independent basis. 
Consequently, the N — 1 dimensional space of vectors orthogonal to W m is spanned by 
the rows of the Gram matrix with the m-th row excepted. The vector (8.4) is orthogonal 
to W m if and only if it lies in the span of the vectors 

(8.5) [z?,...,z%] ne{0,l...N-l}-{m}. 

Thus the condition of Lemma [3] holds if and only if the determinant of the NxN matrix 
whose first (N — 1) rows are the vectors (8.5) and whose last row is (8.4) vanishes as 
stated in the theorem. 

If the weights and the grid points are complex, the same argument can be repeated 
after replacing the weights by their complex conjugates in the definition of W m . □ 

Theorem [4] gives determinantal conditions for boosted order of accuracy. We will cast 
those conditions into a more tractable algebraic form. The following theorem gives the 
template for the algebraic form into which the conditions of Theorem [4] will be cast. 



Theorem 5. If n\,n 2 , 
be factorized as 



,um are distinct positive integers, the determinant (8.3) can 
(zj - Zi) S (zi, ...z N ), 



n 

\<i<j<N 

where S{z\, . . . zji) is a symmetric polynomial that is unchanged when z\,...,zn are 
permuted. All the coefficients of S are integers. 

Proof. We will work over Q, the field of rational numbers. We can think of the deter- 
minant (8.3) as a polynomial in z^ with coefficients in the field Q (z±, . . . , Zjv-i). Since 
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the determinant (8.3) vanishes, if z^ is equal to any one of z±, . . . , Zn-i, we have that 
the determinant can be factorized as 

{z N — Z\) {z N — z 2 ) ■ ■ ■ (zn — zn-i) f 

where / is an element of the field Q (zx, . . . , zjv-i). By Gauss's lemma, / should in fact 
be an element of Z [zi, . . . , Zjv_i], the ring °f polynomials in zx, ■ ■ ■ , Zjv-i with integer 
coefficients (for Gauss's lemma, see Section 2.16 of [T3] and in particular the corollary at 
the end of that section). Now / can be considered as a polynomial in zn-x and factorized 
similarly, and so on, until we get a factorization of the form shown in the theorem. 
To prove that S is symmetric, consider a transposition that switches z v and z q . The 



determinant (8.3) changes sign by a familiar property of determinants. The product of 
all pairwise differences Zj — Z{ also changes sign as may be easily verified or as may be 
deduced by noting that the product is the Gram or Vandermonde determinant. Therefore 
S is unchanged by transpositions and is a symmetric function. □ 

For the determinants that arise as conditions for boosted order of accuracy in Theo- 
rem |4| we describe a method to compute the symmetric polynomial S explicitly. The 
symmetric polynomials that arise in Theorem[5]may well have a connection to symmetric 
function theory. 

To begin with, let us consider the Gram determinant 

(8.6) det {z x , ...,z N , z N+1 ] 0, . . . , N - 1, N) . 
This determinant is equal to 

N 

(8.7) J] ( Zj - Zi) x J] (z N+1 - z k ) . 

l<i<j<N k=l 



See [61 p. 25]. By expanding (8.6) using the entries of the last column (each of these 
entries is a power of Zjv+i), we deduce that the coefficient of z^ +1 in the expansion of 
( |8.6[ ) is equal to 

(8.8) (-i)N+m det ( Zij ^ [o, ... TV — 1, iV] — m). 

This determinant is the minor that corresponds to the entry z^ +1 in the expansion of 



(8.6). By inspecting (8.7), we deduce that the coefficient of z% +1 in that expression is 
equal to 

(8.9) J] - z i) x (-l) N ' m S N - m , 



l<i<j<N 

where 



s P = E 



z il ■ ■ ■ z i p ■ 



l<h<--<i p <N 

Thus S p denotes the sum of all possible terms obtained by multiplying p of the grid 
points Zx, ... z at. For future use, we introduce the notation for the sum of all possible 
terms obtained by multiplying p of the numbers zx, ■ ■ ■ , Zjy, zn+i- 



Theorem 6. The finite difference formula (8.1) with distinct grid points Zf. and weights 



Wk ; m that satisfy (8.2) has an order of accuracy that is boosted by 1 if and only if Sn- 
0. ' 
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Proof. The condition for a boost of 1 is obtained by setting (3 = 1 in Theorem |4| By 
equating (8.8) with ( |8.9 ), we get 

(8.10) det(z 1 ,...,z N ;[0,...N- 1,N] -m) = J[ ( Zj - z { ) x S N - m 

Since the grid points are distinct, the determinant is zero if and only if SN-m = 0. □ 

The corollary that follows covers all the popular cases that have boosted order of 
accuracy. 

Corollary 7. If the grid points Z\, . . . ,zn are symmetric about (in other words z is a 
grid point if and only if —z is a grid point) and N — m is odd, the order of accuracy is 
boosted by 1. 

Although we have restricted m to be in the set {1, 2, . . . , iV — 1}, Theorems [4] and [6] 
hold for the case m = as well. The case m = of (8.1) corresponds to interpolation. 
According to Theorem [6j the interpolation has boosted order of accuracy if and only if 
Sn = or one of the grid points is zero. Of course, the interpolant at zero is exact if 
zero is one of the grid points. We do not consider the case m = any further. 

To derive an algebraic condition for the order of accuracy to be boosted by 2, we 
apply the identity (8.10) with grid points z 1 , . . . , z N , z N+1 and rewrite it as follows. 

det(zi, ...,z N , z N+1 ; [0,...,N-l,N,N+l]-m) = 

N 

II ( Z 3 - Z i) X S N-m+l X - 
l<i<j<N k=l 

We equate the coefficients of z^ +1 to deduce that 
(8.11) 

det(zi,. . . ,z N ;[0, ...N- 1,JV + 1] - m) = Yi (zj - z { ) x (S 1 S N ^ m ~ S N ^ m+1 ) . 

l<i<j<N 

To obtain this identity, we assumed m > 1 and used S N _ m+1 = SV-m+i + "Ziv+iSV-m- 



Lemma 8. The order of accuracy of the finite difference formula (8.1) is boosted by 2 
if and only if SV-m = and SN-m+i = 0. 

Proof. We already have the condition S^-m = for the order of accuracy to be boosted 
by 1. By Theorem [4], the order of accuracy is boosted by 2 if and only if the determinant 
of (8.11) is zero as well. Since SV-m = 0, that is equivalent SV-m+i = 0. □ 



Theorem 9. The order of accuracy of the finite difference formula (8.1) for the m-th 
derivative can never be boosted by more than 1 as long as the grid points are real. Here 
m > 1. 

Proof. By the preceding lemma, the grid points z±, . . . , z^ must satisfy S^-m = and 
Sjv-m+i = for the order of accuracy to be boosted by more than 1. First we consider 
m = 1 and show that Sjv-i — »SW — is impossible. Since iV — 1 > m = 1, we must 
have at least two grid points. Since Sn = 0, at least one grid point must be 0. Since the 
grid points are distinct, no other grid point is zero and iSjy-i ^ 0. 
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If m > 2, let r = JV — m. Then r > 1. To show that S r = S r +i = is impossible, 
denote the elementary symmetric function formed by adding all possible products of r 
numbers out of z±, . . . , Z/v-i by s r . Then 



Si- 



Si- 



S r + Zj^S r -l 
S r+ i + Z N S r 




0. 



Sr— ■ 



Here so is taken to be 1 as usual. Eliminating z^, we get s 2 T 

Newton's inequality (see Theorem 144 on page 104 of [11J) is applied after noting that 
there are at least two numbers in the sequence z±, . . . , ztv_i and that the numbers are 
all distinct. Newton's inequality requires the numbers to be real. We get 

o2 



(V)' 



It is impossible to have s^ 



> 

> 
> 
-is, 



s r+l 



(r-l) (r+l) 

/N -r r + l 
V r N 

Sr— 1 Sr+l 



r — 1 



S r -lS r +l 



-i or S r = S r 



+i 



or S 



N-m 



s 



N—m+1 



0. 



□ 



If the grid points are complex, it may be possible to boost the order of accuracy by 
more than 1. One may obtain formulas for the sequence of determinants with = 1, . . . , b 
in Theore m [I} We have already covered the case with f3 = 1 in (8.10) and the case with 
(3 = 2 in (8.11). To illustrate the general procedure, we show how to get a formula for 
the determinant of Theorem [1] with (3 = 3. We write down the identity (8.11 ) using the 
grid points z\, . . . , zn, zn+i and replace N by N + 1. 



det(zi, ...,z N , z N+1 ; [0, . . . JV - 1, N, N + 2] - m) = 

N 

II ( Z J ~ Z i) X ( S l S N-m+l - S N-m+2) X II ( Z N+1 ~ Z k ). 
l<i<j<N fc=l 

We USe Si = Si + Z N+i , S^_ m+l = S^-m+l + ^Af+l^AT-m, & n d S^f_ m+2 = SN-m+2 + 

ZN+iSjsf-m+i, and equate coefficients of z$ +1 to get 
det(*i, . . . , z N ; [0, . . . JV - 1, JV + 2] - m) = 

IT ( Z 3 ~ Zi ) X — SN-m+lSl + S^-mS^ ~ SjV-mS^) • 

l<i<j<N 

This is the determinant with f3 = 3 in Theorem |4| It gets cumbersome to go on like 
this. However, we notice that the condition for the determinants with /3 = 1,2,3 to be 
zero is SV-m = Sjv-m+i = = 0- Here a simple pattern is evident. 

To prove this pattern, we assume that the determinant of Theorem [4] with /3 = r is of 
the form given by Theorem [5] with 

S = SV-m+r-i + more terms 
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Algorithm 3 Order of Accuracy and Error Constant 



Input: Grid points Z\, . . . , Zn all of which are real. 
Input: Order of derivative m with 1 < m < N — 1. 

Input: Weights it+m, W2, m , • • • , WN,m in the finite difference formula for /^(0). 
Comment: Wk, m are computed using Algorithm 2. 
Input: Tolerance r 



5 



N-m 



T, 



N-m 



Xa<H<- 



,<JV ^ti ■ ■ ■ z 



IN- 



l<ii<—<ijv- 



x<JV 



if 



< rT/v- m then 
— m + 1. 



£>N—m 

r = N 
else 

r = iV — m 
end if 

C = £fc=l w k,mZk 

Leading error term of (1.2): C 



r+m 



/(r+m)(Q) 
(r+m)! 



where each term other than the first has a factor that is one of SV-m, • • • , SiV-m+r-2- 
We pass to the case (3 = r + 1 using the grid points Z\, . . . , z N , z N+ i as illustrated above. 
Then it is easy to see that the form of S for (3 = r + 1 is 

S = S^-m+r + more terms 

where each term other than the first has a factor that is one of S^v-m, • • • > "Sjv-m+r-i- If 
the determinants with = 1, . . . , r in Theorem [4] are zero, the additional condition that 
must be satisfied by the grid points for the determinant with ft = r + 1 to be zero is 

SN-m+r = 0. 



Theorem 10. The order of accuracy of the finite difference formula (8.1) for the m-th 
derivative is boosted by b if and only if «SV-m = Siv-m+i = ■ ■ ■ = SV-m+b-i = 0. Even 
with complex grid points, the order of accuracy can never be boosted by more than m. 

Proof. The first part of the theorem was proved by the calculations that preceded its 
statement. To prove the second part, suppose that the order of accuracy is boosted 
by m + 1. Then we must have Sn = which means at least one of the grid points 
is zero. Since no other grid point can be zero, we must have SV-i 7^ 0, which is a 
contradiction. □ 

Algorithm [3] uses the results of this section to determine the order of accuracy and the 
leading error term in the case of real grid points. We suspect that the error is exactly 
equal to ^ +rj ^ h T for some point ( in an interval that includes and all the grid 
points. 

9. Conclusion 



Algorithm |2j uses partial products of binomials of the type Ilj = i(- 2; — z j) to compute 
finite different weights with good accuracy. A C++ implementation of the method will 
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be posted on the internet. This method lends itself to the efficient computation of 
spectral differentiation matrices as described in Section 5. 

Many finite difference formulas, such as the centered difference formulas for the first 
and second derivatives, have an order of accuracy which is higher than the typical by 1. 
In Section 8, we proved theorems which characterize superconvergence or boosted order 
of accuracy of finite difference formulas completely. 

10. Acknowledgements 

The authors thank John Boyd, Nick Trefethen and Oleg Zikanov for useful discussions. 

References 

[1] R. Baltensperger and M.R. Trammer. Spectral differencing with a twist. SIAM Journal on Scientific 

Computing, 24(5): 1465-1487, 2003. 
[2] J. P. Berrut and L.N. Trefethen. Barycentric Lagrange interpolation. SIAM Review, 46(3) :50 1-517, 

2004. 

[3] F. Bornemann. Accuracy and stability of computing high-order derivatives of analytic functions by 
Cauchy integrals. Foundations of Computational Mathematics, 11(1):1— 63, 2011. 

[4] D. Calvetti and L. Rcichel. On the evaluation of polynomial coefficients. Numerical Algorithms, 
33(1):153-161, 2003. 

[5] R.M. Corless and S.M. Watt. Bernstein bases are optimal, but, sometimes, Lagrange bases are 

better. In Proceedings of SYNASC, Timisoara, pages 141-153. MIRTON Press, 2004. 
[6] P.J. Davis. Interpolation and Approximation. Dover Publications, 1975. 

[7] W.S. Don and A. Solomonoff. Accuracy and speed in computing the Chebyshev collocation deriv- 
ative. SIAM Journal on Scientific Computing, 16:1253, 1995. 
[8] B. Fornberg. Generation of finite difference formulas on arbitrarily spaced grids. Mathematics of 

Computation, 51(184):699-706, 1988. 
[9] B. Fornberg. A practical guide to pseudospectral methods. Cambridge Univ Pr, 1998. 

[10] B. Fornberg. Calculation of weights in finite difference formulas. SIAM review, 40(3):685-691, 1998. 

[11] G.H. Hardy, J.E. Littlewood, and G. Polya. Inequalities. Cambridge Univ Pr, 1988. 

[12] N.J. Higham. Accuracy and stability of numerical algorithms. Society for Industrial and Applied 
Mathematics, 2nd edition, 2002. 

[13] N. Jacobson. Basic Algebra, volume I. Freeman, 2nd edition, 1985. 

[14] B. Sadiq and D. Viswanath. Barycentric Hcrmitc interpolation. Arxiv preprint, 2011. 

[15] C. Schneider and W. Werner. Some new aspects of rational interpolation. Mathematics of Compu- 
tation, 47(175) :285-299, 1986. 

[16] D. Viswanath and L.N. Trefethen. Condition numbers of random triangular matrices. SIAM Journal 
on Matrix Analysis and Applications, 19:564-581, 1998. 

[17] JAC Weideman and SC Reddy. A MATLAB differentiation matrix suite. ACM Transactions on 
Mathematical Software, 26(4):465-519, 2000. 

[18] B.D. Welfert. Generation of pseudospectral differentiation matrices I. SIAM Journal on Numerical 
Analysis, 34(4): 1640-1657, 1997. 
E-mail address: bsadiq@umich.edu and divakar@umich.edu 



